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We study the ground state properties of the interacting spinless fermions in the p^^-orbital bands 
in the two dimensional honeycomb optical lattice, which exhibit different novel features from those 
in the p z -orbital system of graphene. In addition to two dispersive bands with Dirac cones, the tight- 
binding band structure exhibits another two completely flat bands over the entire Brillouin zone. 
With the realistic sinusoidal optical potential, the flat bands acquire a finite but much smaller band 
width compared to the dispersive bands. The band flatness dramatically enhanced interaction effects 
giving rise to various charge and bond ordered states at commensurate fillings of n = | (i = 1 ~ 6). 
At n = |, the many-body ground states can be exactly solved as the close packed hexagon states 
which can be stabilized even in the weakly interacting regime. The dimerization of bonding strength 
occurs at both n = | and |, and the latter case is accompanied with the charge density wave of 
holes. The trimerization of bonding strength and charge inhomogeneity appear at n = |, |. These 
crystalline orders exhibit themselves in the noise correlations of the time of flight spectra. 

PACS numbers: 03.75.Ss,03.75.Nt, 05.50.+q, 73.43.Nq 



I. INTRODUCTION 

There has been tremendous progress during the past 
decade in the cold atom physics. In the early days, Bose- 
Einstein condensation was first realized in magnetic traps 
by using dilute alkali atoms [l|, 0] , where interaction ef- 
fects are weak. Later on, important achievements have 
been made to realize strongly correlated systems by using 
optical lattices. The major advantage of optical lattices is 
the excellent controllability of interaction strength. For 
example, the superfluid to Mott insulator transition of 
bosons has been experimentally observed Recently, 
cold atom physics in optical lattices is merging with con- 
densed matter physics, which provides a wonderful op- 
portunity to explore new states of matter. 

An important aspect of strongly correlated systems 
is orbital physics, which studies an additional degree 
of freedom independent of charge and spin. In many 
transition metal oxides, the d-orbitals are partially filled, 
which enables the orbital degree of freedom active. Or- 
bital physics is characterized by orbital degeneracy and 
spatial anisotropy of orbital orientation. The inter- 
play between orbital, spin and charge degrees of free- 
dom gives rise to many interesting phenomena such as 
metal-insulator transitions, superconductivity, and colos- 
sal magneto-resistance 0, 0, Q ■ 

Orbital degrees of freedom also exist in optical lat- 
tices. Although most of current research of cold fermions 
and bosons focuses on the lowest s-orbital bands, large 
progress has been made in high orbital bands. An im- 
portant advantage of optical lattices is the rigidity of 
lattices. They arc free of the Jahn- Teller type lattice 
distortion which often occurs in transition metal ox- 
ides and quench the orbital degrees of freedom. Orbital 
physics in optical lattices exhibits new features which 
are not usually realized in solid state systems. Recently, 
the properties of bosons in the first excited p-orbital 



bands have been attracting a great deal of attention 
0, 0, H M, HI E GS El- Scarola et al. proposed to 
realize the supersolid state by using bosons in the high 
orbitals to generate the next-nearest neighbor interac- 
tion 0]. Isacsson et al. investigated the sub-extensive 
Zi symmetry of the p-orbital bosons in the square lat- 
tice and its consequential nematic superfluidity [f|. Liu 
and Wu @, and Kuklov [l(| studied the antiferromag- 
netic ordering of the on-site orbital angular momentum 
moment. It was also proposed in Ref. 0] to enhance the 
life-time of p-orbital bosons by using a Bose-Fermi mix- 
ture to reduce the available phase space of decay process 
of bosons. Wu et al. [ll[ further investigated the su- 
perfiuid and Mott insulating states of p-orbital bosons in 
the frustrated triangular lattice, and found a novel stripe 
phase of orbital angular momcntums. Xu et al. studied 
a model of bond algebraic liquid phase [l2j and phase 
transitions in anisotropic xy-models fl3| in the context 
in the p-orbital boson systems. 

On the experimental side, the progress of orbital 
physics with cold atoms has also been truly exciting, 
which opens up the new opportunity to stud y o rbital 
physics. Browaeys et al. [15| and Kohl et al. [16| have 
demonstrated the population of high orbital bands with 
both bosons and fermions. Furthermore, Sebby-strabley 
et al. fl7| have successfully pumped bosons into the ex- 
cited bands in the double-well lattice. More recently, an 
exciting progress has been made by Mueller et al. [HI 
to realize the meta-stable p-orbital boson systems by us- 
ing the stimulated Raman transition to pump bosons to 
high orbital bands. The spatially anisotropic phase co- 
herence pattern has been observed in the time of flight 
experiments. This opens up a new experimental direction 
to investigate novel condensate of bosons in the excited 
p- bands. 

On the other hand, fermions in the p-orbital bands 
also possess interesting behaviors [lj| [2(| Hi], [22]. Re- 
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cently, Wu et al. [19( studied the flat band structure in 
the Pz.y-orbital physics in the honeycomb lattice. Com- 
pared to the p 2 -orbital system of graphene, which has 
been attracting tremendous attention since the discov- 
ery of quantum Hall effect therein [23], [U [25| , the p x ,y- 
orbital honeycomb systems exhibit new and even richer 
physics. In graphene, the active bands near the Fermi 
energy are 'V'-type, composed of the p z -orbital directly 
normal to the graphene plane, thus graphene is not a 
good system to investigate orbital physics. In contrast, 
it is the other two p-orbitals {p x ,y) that lie in-plane and 
exhibit both orbital degeneracy and spatial anisotropy, 
giving rise to the interesting flat band physics (l9j |. In 
solid state systems of graphene and MgB 2 , p X)J ,-orbitals 
hybridize with the s-orbital, resulting the er-bonding (the 
sp 2 hybridization) band. This er-band is fully filled and 
inert in graphene, but is partly filled and contributes to 
the two-band superconductivity in MgB2 [26[ . Due to the 
large s-orbital component, the essential feature of orbital 
physics, orbital anisotropy, is not prominent in these two 
systems. In contrast, the p^y-orbital bands in optical lat- 
tices are well separated from the s-band with negligible 
hybridization, providing an unique opportunity to study 
the pure p^-orbital physics in the honeycomb lattice. 
This research will provide us another perspective in the 
honeycomb lattice and is complementary to the recent 
research focus on the single band system of graphene. 
Other works of the p-orbital fermions include the investi- 
gation of orbital exchange physics in the Mott-insulating 
states finding various orbital ordering and frustration be- 
havior 0, [22[, and the study of the possibility to en- 
hance the antiferromagnetic ordering of fermions in the 
p-orbital of 3D cubic lattices pnj . 

Interaction effects in the p^^-orbital honeycomb opti- 
cal lattices can be much stronger that those in the p z - 
orbital graphene systems. In real graphene the dimen- 
sionless coupling constant r s = e 2 /(eKv) has a maximum 
value of 2.3 in vacuum (and r s < 1 for the current avail- 
able graphene samples on SiC-2 or SiC substrates), tak- 
ing v = 10 6 cm/sec. Thus graphene is very far from 
the r s = 39 regime needed for Wigner crystallization 
[2^ . Much of graphene interaction physics is described 
by perturbative weak-coupling renormalizations of the 
quasiparticle spectral function, as shown both theoret- 
ically and experimentally (28l. I2DI |30| . Furthermore, real 
graphene physics is complicated by electron-phonon in- 
teractions 13 1| . In contrast, in the p^^-orbital honey- 



comb lattices systems, the flat band quenches the kinetic 
energy, and thus interaction physics is non-perturbative 
and generic, leading to qualitatively new orbital physics 
phenomena, e.g. Wigner-Mott physics, can show up eas- 

ily [H- 

This paper works as an expanded version of a previous 
publication of Ref. [lj| , with new results and all the the- 
oretical details of the behavior of spinless fermions in the 
Pa;,i/-orbital bands in the honeycomb lattice. The current 
work is motivated by considerations of using the optical 
lattices to go beyond what can be achieved in solid state 



systems, i.e. obtain exotic strongly correlated orbital 
quantum phases which have not yet been studied in con- 
densed matter physics. The paper is organized as follows. 
In Sect. HH we analyze the band structures in both the 
simplified tight-binding model and the realistic optical 
potential constructed from three co-planar laser beams 
[32l | . The band structures contain both Dirac cones in two 
dispersive bands, and other two nearly flat bands over the 
entire Brillouin zone (BZ) whose flatness becomes exact 
if the 7r-bonding is neglected. Special attention is paid 
for the orbital configurations of the localized Wannier- 
like eigen-functions in the flat bands, and also at the 
Dirac points. In Sect. IIII1 the interacting Hamiltonian 
is introduced and methods of enhancing the Hubbard- 
like on-site interaction are proposed. In Sect. IIV1 the 
interaction effect in the partially-filled flat band is dis- 
cussed. The situation is somewhat analogous to that in 
the fractional quantum Hall effect of electrons in the low- 
est Landau level. When the flat band is partially-filled, 
the effects of interactions are entirely non-perturbative. 
We obtain the exact many-body plaquette Wigner crys- 
tal state at filling (n) — i, which is the close packed 
hexagon state and is stable even in the weak interaction 
regime. In Sect. [V] we present various charge and bond 
ordered states, including dimerized and trimerized states 
at higher commensurate fillings in the strong interaction 
regime. In Sect. IVI1 the noise correlation in the time of 
flight experiments are discussed. Conclusions and out- 
look for the future research are discussed in Sect. IVIII 



II. P:c ^-ORBITAL BAND STRUCTURE IN THE 
HONEYCOMB LATTICE 



In this section, we will give a detailed analysis to the 
Pa;,j/-orbital band structure in the 2D honeycomb lattice 
which is featured by the interesting properties of both flat 
bands and Dirac cones. We will first discuss the experi- 
mental construction of such a lattice and then solve the 
band structure by using both the simplified tight-binding 
model and the realistic sinusoidal optical potential. 



A. Construction of the optical honeycomb lattice 

The honeycomb optical lattice was realized experimen- 
tally by using three laser beams with co-planar propa- 
gating wavevectors cft(i = 1 ~ 3) quite some time ago 
[32j . The magnitudes of these wavevectors are the same 
and their directions form the angle of 120° with each 
other. Assuming the polarization of the electric fields of 
the three beams are all along the z-direction, the optical 
potential distribution can be expressed as 



V(r) = V ^2 C0S (P* ' 



(1) 



where pi = q 2 - q 3 , p 2 = q 3 - qi, and p 3 = q\ - q 2 . In 
the case of blue detuning, Vq is positive and the potential 



FIG. 1: A) The contour plot of the optical potential of the 
2D honeycomb lattice described by Eq. Q] B) The optical 
potential distribution around the potential maximum in one 
unit cell. 
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FIG. 2: A) The two sublattice structure (A and B) of the 
honeycomb lattice. B) The hexagon Brillouin zone with edge 
length 4?r/(3v / 3a). 



minima form a hexagonal lattice as depicted in Fig. [T] 
In contrast, the red detuning laser beams generate a 2D 
triangular lattice. In both cases, the lattice is topolog- 
ically stable against the phase drift of the laser beams, 
which only causes a overall shift but not the distortion of 
the lattice. Fig. Q]B depicts the potential distribution in 
one unit cell of the honeycomb lattice, where a potential 
maximum locates in the center and six potential min- 
ima sit around. Without loss of any generality, we take 
P1.2 = p(±^e x + \e y ) and p 3 = -p e y where p = |^ 
and a is the distance between the nearest neighbour site 
in the honeycomb lattice. We define the recoil energy in 
such a lattice system as E r — where m is the mass 
of the atom. 



B. The tight-binding model 

The optical potential around the center of each site 
is approximately an anisotropic harmonic well. We as- 
sume that the vibration frequencies along the x, y and 
z-directions satisfy u z S> uj x = uj y = u> X y, and thus the 
energy of the p z -orbital is much higher than that of the 
p x ,y-orbit&\ bands. When the lowest s-band is fully filled 
and thus inert, the active orbital bands will be of the 
p x ,y Due to the spatial orientation of the p-orbitals, the 
hopping processes in the p-orbitals can be classified into 
the a and 7r-type bondings, respectively. The former de- 



FIG. 3: Dispersion of the two-lowest p^-orbital bands E\^,- 
The band E\ is completely flat, while E2 exhibits Dirac points 
at K\,i = (zfc , 0). The other two bands are symmetric 
with respect to E = 0. From Wu et al 



scribes the hopping between p-orbitals on neighbouring 
sites with the orientation along the bond direction, while 
the latter describes the hopping between p-orbitals per- 
pendicular to the bond direction. In other words, the <j- 
bonding is of the "head to tail" type, while the 7r-bonding 
is of the "shoulder by shoulder" type. Typically, the am- 
plitude of the 7r-bonding is much smaller than that of the 
cr-bonding because of the strong orientational anisotropy. 

The structure of the honeycomb lattice is depicted in 
Fig. [2]A. Each unit cell in the honeycomb lattice contains 
two sites depicted as A and B. We define three unit 
vectors from site A to its three neighbouring sites B as 



ei,2 = ± ~^ e x 



1 



and their differences 6, = -ke 



2 v ' 

V3 



v3e a , &1.2 = -— e x ± -e y . 



(2) 



(3) 



The projections of the p-orbitals along the &i,2,3 direc- 
tions are defined as 



V3 1 

Pl,2 = ±— Px + 7,Py> 



P3 



-Py 



(4) 



Only two of them are linearly independent. In the real- 
istic optical potential depicted in Fig. Q] A, the poten- 
tial distribution inside each optical site is only approxi- 
mately isotropic in the sy-plane. Away from the center, 
the potential exhibits a 3-fold rotational anisotropy. The 
point group symmetry respect to the center of each site 
is reduced into Czv including the 3-fold rotation and re- 
flection. Nevertheless, as required by this symmetry, p XtV 
remain degenerate and each of pi,2.3 defined above is still 
the on-site eigenstates with the orientation along the cor- 
responding bond direction. But they are no longer purely 
parity odd due to the breaking of the inversion symme- 
try respect to the center of each optical site. (The overall 
inversion symmetry respect to the center of the each po- 
tential maximum is still preserved, but this involves the 
transformation among different sites.) 
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The (T-bonding part in the kinetic energy reads 



Hn 



+aei ,1 



h.c.} — (i n?, (5) 

regies 



where the summation over r in the first term is only on 
the A sublattice, a is the nearest neighbor distance, and 
np = np lX + np, v is the total particle number in both p x 
and p y orbitals at the site r. t» is positive due to the dom- 
inant odd parity component of the p-orbitals and is set 
to 1 below. Eq. [5] neglects the much smaller 7r-bonding 
t± terms which in principle exist, and their effects will 
be discussed in Sect. Ill El 

Next we discuss the spectrum of the tight-binding 
Hamiltonian Eq. [5] In momentum space, we define a 
four-component spinor as 



ip(k) = (pAx(^),PAy(fc),ps a; (fc),p By (fc)) 5 



(6) 



where each component is the Fourier transform of the 
Px,y-OTbit in site A or B. Then Eq. [5] becomes 



H 



*|lI>£(fe){jMfc)-M<W}^(fe), (?) 



where the matrix kernel H a p(k) takes the structure as 



/ f(e'*- ffl +e ife ' ?2 ) 



V 




h.c. 








— (e ife?1 — e ik ' 32 ) 






) 



Its spectrum is symmetric respective to zero because 
the sign of the in term can be flipped by changing the sign 
of the p^y-orbitals in one sublattice but not the other. 
The dispersion relations of the four bands read 



3 *|| / 

#1,4 = T 2*11 ' #2,3 = /3- 



2 ^ cos • 6j (8) 



as shown in Fig. [3l Interestingly, the band structure 
exhibits two flat bands over the entire 2D Brillouin 
zone. The corresponding eigenvectors can be found ana- 
lytically as 



■0i,4(fc) 



1 



N (k) 



{^[/ 2 * 3 (fc)-/ 3 *i(% -/r 2 (fc), 



± 



k/23(fe)-/3l(fc)], T/l2(fc)} T , 



where fij 
reads 



V3 

ik- Ej 



(9) 



e ik-ej anc j ^- ne norma li za tion factor 



N (k) = -(3-^cosfc-^). 



(10) 



On the other hand, the E2.3 bands are dispersive ex- 
hibiting the Dirac cone structure, whose band width is 




FIG. 4: The Wannier-like localized eigenstate for the lowest 
band. The orbital configuration at each site is oriented along a 
direction tangential to the closed loop on which the particle is 
delocalized. The absence of the 7r-hopping and the destructive 
interference together ensure such a state as an eigenstate. 



determined by fn . We construct a new set of basis which 
are orthogonal to ipi^k) and span the subspace for the 
#2,3 bands as 

4>{k) = \Jj^{fa&)> ^(/23(fc)-/3l(fc)), 0, 0}, 

4>'{k) = ^{0, 0, Am ^(/2 3 (£)-/ 3 *i (£))}• 

(11) 

Then the Hamiltonian becomes the same as in graphene 



#2 3 (fc) = 



^ — ik-ei 



2 VE, 



(12) 



Two Dirac cones appear at K 12 = (± 3 ^ ,0). The 
eigenvectors of the bands #2,3 (fc) read 

^ 3 (k) = ^={cb(k)±e l9 ^'(k)}, (13) 

with the angle of B-^ 

fe = arg£V^). (14) 

i 

C. The localized eigenstates 

The complete flatness of the £1,4 bands means that 
these eigenstates can be represented as linear superposi- 
tion of a set of degenerate localized states. The construc- 
tion of these localized state is depicted in Fig. SJ For 
each hexagon plaquette denoted by its center position R, 
there exists one such eigenstate for the bottom band E\ 



6 

E 



(-Y 1 j cos^lp^) - sm6j\p jty )j, (15) 
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where j is the site index and 9j = {j — The lo- 

calized eigenstates of the E4 band can be obtained by 
flipping the signs of the p-orbits on sites 2, 4 and 6 and 
keep those on sites 1, 3 and 5 unchanged. The p-orbital 
configuration on each site is perpendicular to the links 
external to the hexagonal loop, thus the a-bonding for- 
bids the particle to directly hop outside through these 
links. Furthermore, the amplitudes for the particle to 
hop to the p-orbital in the radial direction from the neigh- 
bouring sites vanish due to the destructive interference 
as shown in Fig. 3J The particle is trapped in the pla- 
quette without "leaking" to outside, and thus \ipg) is the 
eigenstate with the energy of E\. The states are 
all linearly-independent apart from one overall constraint 
I^r) ~ under periodic boundary conditions. The 
localized states on two neighbouring edge-sharing plaque- 
ttes are not orthogonal to each other. 

The Bloch wave states in the flat band E\ are con- 
structed as 




The doubly degenerate eigenstate at k = (0, 0) can not 
be constructed from the above plaquette states. They 

are l^fe=(0,0))l.2 = EreA \Px(v),r) ~ J2feB \Px{y),r)- 




B) % 



FIG. 5: The orbital configurations of eigenstates at Ki, which 
are of p x ± ip y type. The phase of each lobe is presented 
(w = e z2 r ). A) Vi(^i) (Ei = B) MKi) with a k = 

(Ba = 0). 



D. Orbital configuration at k = (0,0) and Ki t 2 

The major difference between the physics of p x ,y- 
orbital bands and that of graphene is the orbital degree of 
freedom. The orbital configuration for each band varies 
as lattice momentum k changes in the Brillouin zone. 
Around the center of the Brillouin zone k = (0,0), the 
Hamiltonian can be expanded as 



K\ , the Hamiltonian can be expanded as 
H(k) = 



3 3 
-Afc^Ti <g> I + -Ak v r 2 <g> I 

4 4 



3 3 3 

(- + -Ak x ) n <g> ct 3 - -Ak y t 2 <8> cr 3 

4 8 8 

3 3 3 

-Aky n <g> o-i - (- - -Ak x )r 2 ® ax, (19) 



H = -ti ® I - -k y r 2 ® ct 3 - -k x T 2 O cri, (17) 

where Pauli matrices CTi.2,3 describe the p^.y-orbital de- 
grees of freedom, and 71,2,3 describes the sublattices A, B 
degrees of freedom. The eigenvectors of ^1,2,3,4 can be 
approximated as 

■0i,4(fc) = ^ { - k y ,k x ,±ky,Tkxj, 



i>2,3(k) = -^^k x ,k v ,Tk x ,Tky}- 



(18) 



Thus around k = (0, 0), the orbital configuration is polar- 
like, i.e., a real combination of p x ,y The orbital orien- 
tation in each site is either parallel or perpendicular to 
k. 

Now let us investigate the orbital configurations 
around the vertices of K\ 2 of the Brillouin zone. Around 



where Ak = k — K\ and g±(k) = Ak x ±iAk y . The eigen- 
vectors of the flat bands E\ 4 can be approximated as 



Mk) = i{l + 



ff+O) g+(fc) j 



± n+tM], ±i[i 



} • (20) 



Similarly, the eigenvectors of the dispersive bands E 2 ^ 
are approximated as 

±e-(l-^, -i(l + ^)} T ,(21) 

where is the angle defined in Eq. Q3] Thus the orbital 
configuration at k = K\ on each site is the axial state p x ± 
ip y as depicted in Fig. O This is in contrast to the polar 
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FIG. 6: Dispersion of the two-lowest p^.y-orbital bands E\ t 2 
in the presence of the 7r-bonding t±. The bottom band is no 
longer rigorously flat but acquires a narrow width of t±. 




FIG. 7: Dispersion of the two-lowest p^y-orbital bands E1.2 
in the presence of the on-site potential. 

configuration at k = (0,0). The orbital configuration 
at k = K 2 can be obtained by performing time reversal 
transformation. 

E. 7r-bonding term and other perturbations 

The 7r-bonding term in principle exists in the realistic 
optical lattices. We define the projections of p^-orbitals 
perpendicular to the £1,2,3 directions as 

Pi, 2 = ~\px ± -^-Py P' 3 =Px- (22) 
The 7r-bonding term can be written as 

= ~t X {P'hp'r+ae^+h.c.}. (23) 

reA,i=l~3 

Please note that the hopping integral of the 7t-bonding 
has the opposite sign to that of the tr-bonding. In this 
case, the bottom and top bands £1,4 are no longer rig- 
orously flat but develops a narrow width 3t± as depicted 



in Fig. [6] The E\ and E 2 bands still touch at the center 
of the Brillouin zone. This can be understood from the 
structure of the localized eigenstates in Fig. 0J The Tr- 
bonding term causes the particle leaking off the plaquctte 
and thus correspondingly develops the band width. Nev- 
ertheless we will show in Sect. Ill Fl that in the realistic 
optical potential that such an effect is negligibly small. 

Next we discuss the case that the A and B sites are 
with different on-site potentials. In the graphene-like sys- 
tems, this corresponds to a mass term in the Dirac point. 
In the p^y-orbital systems, such a term can be described 
as 

H = AE^J^np-J^npy (24) 

r£A reB 

The spectrum is depicted in Fig. [7] with the opening of 
a gap of AE in the Dirac points as usual. Interestingly, 
the flat band feature remains unchanged. This can be 
understood in terms of the localized eigenstate picture. 
In this case, the localized eigenstates of the E\ band still 
possess a similar configuration as in Fig. 01 but their 
wavefunctions distribute with different weights on A and 
B sublattices. 



F. Band structure from the continuum optical 
potential 

We numerically calculate the band structure in the re- 
alistic optical potential of Eq. [TJ The band Hamiltonian 
becomes 

H = - ^ + V £ cos® • r). (25) 

i=l~3 

Since this is a non-singular sinusoidal potential, we use 
the plane wave basis to calculate the matrix elements 
(k\H\k') where k' = k ± = 1 ~ 3). For each k 
in the Brillouin zone, we truncate the matrix up to 120 
plane- wave basis, which should be sufficient for the lowest 
several bands. 

The band dispersions along the path from O to K\ , M 
and K 2 are depicted in Fig. [8] The locations of O, ^1,2, 
M in the Brillouin zone are depicted in Fig. [51 The lowest 
two bands are of the s-orbital exhibiting Dirac cones at 
K\ and K 2 . The next four are of the p^^-orbitals. The 
band flatness is largely preserved even with the realistic 
optical potential of Eq. [Q In Fig. [5] A with V/E r = 5, 
the bottom one of the four p-orbital bands is nearly flat 
with the width of 7 x 10~ 3 £V which is only 2% of that 
of the second one which is 0.35E 1 ,-. The width of the top 
band is 4 x 10 -2 -E r which is still small but considerably 
larger than that of the bottom one. The third band is 
the widest one with the width of 0.62E r . As we can see, 
the particle-hole symmetry in the tight-binding model is 
no longer kept because of the unavoidable hybridization 
with other bands and long range hoppings. The spectra 
become more symmetric with a strong optical potential 
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FIG. 8: Band structures for the realistic optical potential in 
Eq. Q] along the path from O — + K\ — > M — > A"2 in the 
Brillouin zone with A) V/E r = 5 and B) V/E r = 10. 



as shown in Fig. [8] B (V/E r = 10) in which the tight- 
binding model is a better approximation and long range 
hoppings can be neglected. The widths of the beginning 
six bands (the s and p^-orbital bands) as a function of 
the Vo/E r are depicted in Fig. O 



III. INTERACTIONS IN THE p^-ORBITAL 
SPINLESS FERMIONS 

In the following, we will mainly consider interact- 
ing spinless fcrmions in the p x ^-orbital bands in the 
honeycomb lattices, and leave the research for spinful 
fermions in future publications. The preparation of spin- 
less fermions can be controlled by cooling the system in 
the external Zeeman field. Due to the lack of spin relax- 
ation mechanism in cold atom systems, the system will 
remain in the spin polarized state. The spinless fermions 
have been realized in many experiments. In particular, 
the strongly correlated polarized spinless fermion systems 
have been realized by using the p-wave Feshbach reso- 
nance [33|, H3, [H, [3y|- Therefore, in contrast to solid 
state electronic systems, where spin is almost always an 
important quantum dynamical variable, the cold atom 



FIG. 9: Width of the beginning six bands of the optical po- 
tential Eq. [1] Bands 1 and 2 are the s-orbital bands. And 
bands 3,4, 5, and 6 are the p x<y orbital bands, where those of 
3 and 6 are flat as discussed before. 

fermionic systems created by Feshbach resonance, can 
be prepared as spinless (i.e. spin polarized), and conse- 
quently, our current spinless theory applies to such sys- 
tems without any modifications. Of course, the problem 
of creating a laboratory p x y orbital graphene system in 
cold atomic gases still remains, but given the rapid cur- 
rent experimental developments in fermionic cold atom 
matter, we are optimistic that our prposed system should 
soon be realized in practice. 

Because of the orbital degeneracy, the on-site interac- 
tion for spinless fermions remains Hubbard-like 

Hmt = U y^^n^xn^y, (26) 

r 

where the on-site interaction U is 

U = J dfidftjV^i - r 2 ) ((^ {ri)^ Py (r 2 )) 2 

- (^>l)Vv>l)Vv>^ Py (r 2 )). (27) 

Due to Paul's exclusion principle, the s-wave scattering 
vanishes, and thus the p-wave scattering is the leading 
order contribution which is typically weak for low energy 
particles. The p-band fcrmions have high kinetic energy, 
and thus their p-wave scattering might not be small. To 
enhance U, we can use the p-wave Feshbach resonances 
among spinless fermions (e.g. 40 K [H, UK], 6 Li[35|). Al- 
though we do not want the system staying too close to 
the resonance because of the large atom loss rate there, 
an enhancement of U to the order of the recoil energy 
En while maintaining the stability of the system is still 
reasonable. 

Another possible method is to use atoms with large 
magnetic moments which interact through magnetic 
dipole-dipole interactions as 

V(n - f 2 ) = -^{mi ■ m 2 - 3(toi ■ f)(m 2 • f)}, (28) 
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where r = \r± — f^l and f = (ri — r^jr. The fermionic 
atom of 53 Cr is a good candidate whose magnetic moment 
is mcr = GfJ-B (Bohr magneton). The spin polarization 
can be controlled by an external magnetic field. Below 
we give an estimation of U from the magnetic dipole in- 
teraction. The vibration frequency in each site can be ob- 
tained as u! x .y = \J §Vo£y. The length scale of the peri- 
orbitals (l x . y — yjh/muj x ^ y ) is typically one order smaller 
than a. For example, we estimate that l x , y /a » 0.2 at 
Vo/E r = 5. Assuming strong confinement in the z-axis 
h <C lx,y> the vector fx — r- 2 linking two atoms in p x and 
p y orbits almost lies in the xy-plane. When the fermion 
spin is polarized along the z-axis, the interaction is re- 
pulsive and U can be approximately estimated as 

9 

977 ^ 

U^^(l-3(co S 2 6)), (29) 

where 9 is the angle between f\ — fj and the z-axis. We 

estimate (r) = yj 2l^ y + 1\ and cos# = l z /r, and find 

that U can reach the order of E r . For example, if we 
use the laser wavelength A w 0.8/zm, V/E r = 5 (so that 
l x y k, 0.2a) and l z = 0.2l x y , we arrive at U — 2.2 KHz or 
approximately lOOni'T. Increasing V/E r can further in- 
crease U and suppress in , thus drive the system into even 
stronger correlation regime. U can be adjusted from re- 
pulsive to attractive by tuning the polarization direction 
from perpendicular to parallel to the xy-plane. 



IV. THE WIGNER CRYSTAL STATE AT 
i-FILLING 

6 

In this section, we discuss interaction effects in the 
Px,j/-orbital systems. When the flat band is partially 
filled, interaction effects dominate the physics. In par- 
ticular, a Wigner crystal state is stabilized even with the 
shortest range on-site interaction. We will mainly study 
the spinless fermion system below, and also give a brief 
discussion on the boson systems, but leave the study of 
the spinful fermion systems to a future publication. 

A. Close packed plaquette state 

Due to the complete suppression of the kinetic en- 
ergy in the flat band, the effect of interactions is non- 
perturbative when the flat band is partially filled. Inter- 
estingly, at sufficiently low particle density n < 4, the 
exact many-body ground state can be easily constructed 
as follows. Each individual particle localizes into a pla- 
quette state depicted in Fig. [4] Any arrangement of 
these plaquette states avoiding touching each other is 
the kinetic energy ground state and costs zero interac- 
tion energy. Since the interaction is repulsive, this class 
of states also minimize the interaction energy and thus 
they constitute the many-body ground states. If we fix 




FIG. 10: The configuration of the close-packed Wigner- 
crystal state for both bosons and fermions at n = 1/6. Each 
thickened plaquette has the same configuration as in Fig. 3] 
From Wu et al 

the particle density at n < 1/6, the ground state config- 
urations have large degeneracy corresponding to all the 
possible ways to arrange these hard hexagons. 

Another class of systems exhibiting similar behavior is 
the frustrated magnets near full polarization in a large ex- 
ternal magnetic field. The Holstein-Primakoff magnons, 
which are bosons, have a dispersionless flat band over the 
magnetic Brillouin zone. Interactions among magnons 
result in the magnon crystal state and magnetization 
plateau [37j near the full polarization. However, this 
flat band behavior is difficult to observe because a very 
strong magnetic field to drive the system close to the full 
polarization is required. This means that the Zeeman 
energy reaches the exchange energy J which is typically 
larger the order of meV. Flat band phenomenon also ap- 
pears in systems of "fermion condensation" where strong 
interactions drive an originally dispersive band to flat 
within a finite width around the Fermi energy 38]. This 
has been proposed to explain the Curie's law behavior of 
the magnetic susceptibility in the itinerant heavy fermion 
compound CeCoIn 5 system [o9| . 

The close packed plaquette pattern without overlap- 
ping each other is depicted in Fig. [TO] corresponding to 
the filling of n — |. The completely filled lowest flat 
band corresponds to n = |, thus this close packed pla- 
quette pattern corresponds to i-filling of the flat band. 
This state breaks the lattice translational symmetry and 
is three-fold degenerate. The other two equivalent states 
can be obtained by translating the state in Fig. [TO] along 
x-axis in the right or left direction at one lattice constant. 

B. Stability of the Wigner crystal state 

The above Wigner crystal state is a gapped state. We 
can give a rough estimation for a upper limit of the charge 
gap by constructing a trial wavefunciton for putting an 
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i-1.3- 




range of 



~ 1. The self-consistent equation reads 



FIG. 11: The phase boundary of the incompressible plaquette 
Wigner crystal state of spinless fermions at (n) — |. From 
Wu et a/ [g. 



extra particle in the close packed state in Fig. |TD] In 
the weak interaction case {U/t\\ <C 1), we can put the 

extra atom in the plaquette state located at R which 
is adjacent to three occupied plaquettes -Ri.2,3- Since 
there is already 1/6 atom on average per site, the cost 
of the repulsion is ¥. On the other hand, in the strong 
coupling case (U/t\\ ^S> 1), we put the particle into an 

excited state of the occupied plaquette R\ while fixing the 
orbital configuration on each site. Because fermions are 
spinless, the cost of energy comes from the kinetic part 
with the value of . Thus we obtain the upper limit 
for the charge gap which is determined by interaction at 
small values of U and by kinetic energy in at large values 
of U as 



A<min(i/J,~t||). 



(30) 



(n? ti ) = (fi[(7V,j)] \n?,i\ £}[vV,j)]) (ij = 1 ~ 3), (33) 

where |0[riffj-]) is the mean field ground state with the 
specified configuration of (npj). 

At the mean field level, we found that (np 3 ) is zero 
which means the time reversal symmetry is kept. We 
need to take an enlarged unit cell to allow the spatial 
variation of the order parameters. In order to obtain the 
plaquette order in Fig. [TOI this enlarged unit cell covers 
six sites around a plaquette. We present the range of 
chemical potential [i for the i -Wigner crystal state in 
Fig. [TTJ which corresponds to the excitation gap. The 
charge gap grows roughly linearly with U in the weak 
interaction regime, and saturates at a value comparable 
to t\\ in the strong interaction regime. Both agree with 
the above variational analysis. 

Next we discuss the effect of the 7r-bonding term to 
the g-state. Such a term gives a width determined by 
t± to the originally flat band. Because the g-state is 
gapped, it should remain stable if the band is sufficiently 
narrow. The plaquette state costs the kinetic energy at 
the order of t± per particle while it saves the repulsive 
interaction at the order of ^ . Thus for small values of t± , 
a stability condition of this state can therefore be roughly 
estimated as U > 6t± . We have checked this numerically. 
For example, setting the ij_/f|| = 0.1, we find that the |- 
state survives U > t» . In realistic systems, the ratio of 
t±/tii is much smaller than 0.1 with reasonable values 
of V/E r as shown in Fig. [9j thus the 1/6-state can be 
stabilized at much smaller values of U . 



The above intuitive picture can be make more rigorous 
by performing self-consistent mean field treatment to the 
Hamiltonian described below. 

We decouple Eq. [55] both in the direct and exchange 
channels as 



H, 



mf,int 



= u 



— pt V- 



{n?,i) - i(np,3) 



h.c 



}, (31) 



np,i = 




PpyPry), 




n?,2 = 


2^-PfxPry ~ 


f h.c), 




?V,3 = 


2~l (PrxPry 


- h.c), 


(32) 



where ni,2,3 are the pseudo-spin operators. np,i,n?^ are 
time-reversal invariant, and describe the preferential oc- 
cupation of a "dumbbell-shaped" real p-orbital orienta- 
tion; n^3 is the orbital angular momentum, and is time- 
reversal odd. We perform a self-consistent mean field 
solution to Eq. (3D plus Eq. [5] for the filling level in the 



C. Bosonic Wigner crystal state 

In the above hard hexagon state at n = ^ , particles are 
separated from each other, thus particle statistics do not 
play any role. Such a Wigner crystal state should also 
occur with bosons or Bose-Fermi mixtures with repulsive 
interactions. The on-site interaction for p-band bosons 
reads 



Hint — 



(34) 



where n is the total particle number and L z is the or- 
bital angular momentum fill ]. The p-band bosonic 
systems have been created experimentally [l8j]. The life 
time of the p-band bosons is significantly enhanced when 
the particle density per site is less than one, which can 
be hundreds of times longer than the hopping timer from 
one site to its neighbours. Thus the 1/6-state is also 
experimentally feasible in bosonic systems. 
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FIG. 12: The filling (n) vs. the chemical potential /1 for 
spinless fermions for weak A) and strong B) interactions. Due 
to the particle-hole symmetry, only the part with fi from the 



to U/2 is shown. Only one plateau appears 

ear in B) at 



band bottom — 1 1\ 

in A) at n = |, while a series of plateaus ar. 
n = 1/6, 1/3, 1/2, 2/3, 5/6, 1. From Wu et al 1 



V. CHARGE AND BOND ORDERINGS AT 
COMMENSURATE FILLINGS OF n > ± 

In this section, we investigate the charge and bond or- 
derings at commensurate fillings higher than g by using 
the mean field theory to solve the interacting Hamilto- 
nian self-consistently. We will present the result in both 
weak and strong coupling regimes, but leave the detailed 
investigate of the physics of orbital exchange at n = 1 to 
a future research. In the following calculation, we confine 
ourselves to the unit cell up to 6-sites. 



A. Weak coupling regime 

When the filling n > 1/6, exact solutions are no longer 
available. Again we perform the self-consistent mean field 
solution to the interacting Hamiltonian. In the weak cou- 
pling regime (U/ti\ = 1), we plot the relation of the filling 
(n) vs. [i in Fig. [T^] A. As /j, passes the charge gap, the 




FIG. 13: Bonding strength dimerization can occur at both 
(n) = 1/3 in the weak coupling regime and (n) — 1/2 in 
the strong coupling regime as depicted by the thickened (red) 
bonds. The orbital orentation in the dimer is along the bond 
direction. In the weak coupling case ((n) = 1/3), the the 
thickened bonds correspond to the shared edges of two neigh- 
boring plaquette states in the flat band. In the strong cou- 
pling case ((n) = 1/2), each dimer contains one particle as an 
entangled state of occupied and empty sites. 



system enters a compressible state, (n) increases with /1 
quickly with a finite but large slope. This means that 
particles fill in other states in the flat band. Due the 
preexisting crystalline ordered background, these states 
are no longer exactly flat and develop weak dispersions. 
This corresponds to adding additional particles into 
the background of the i -state as depicted in Fig. [TU] 
(N: the total number of lattice sites). Roughly speak- 
ing, these new particles also go into the localized plaque- 
tte states. When (n) w |, we see a significant reduction 
of density of states compared to those of the flat bands, 
which still has a finite density of state attributed to the 
filling of the dispersive band. 



Let us look at the quasi-plateau at 



i Although 



these newly occupied plaquettes can be arranged to avoid 
each other as we did before, they unavoidably will touch 
the preoccupied ones. As a result, for each occupied pla- 
quette state, three of its six neighbours are occupied al- 
ternatively. The orbital configuration in such a state is 
as depicted in Fig. Q21 for each bond shared by two 
occupied plaquettes, the p-orbital orientation is paral- 
lel to the bond direction as a compromise between two 
neighbouring plaquettes. The bonding strength exhibits 
a dimerized pattern. The ratio between the weakened 
and strengthened bonds is approximately 0.44. Com- 
pared to the gapped dimerized phase discussed below in 
Sect. IV Bl this is a relatively weakly dimerized phase. 

At (n) > 1/2, all the flat band states are completely 
filled. In the weak coupling regime, interaction effects are 
no longer important and crystalline orders vanish. Near 
(n) = 1, the density of states becomes linear with energy 
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FIG. 15: The orbital configuration at the filling (n) =5/6 
exhibits the dimerized state of holes in the strong coupling 
regime. Holes mainly distribute on the position of sites 1 and 
4 in each unit cell with a large bonding strength. 

as controlled by the Dirac cones. Recently, it has been 
proposed to use the s-band in the honeycomb optical 
lattice to simulate the Dirac cone physics 40]. The p x ,y- 
band Dirac cones described above are also good for this 
purpose and have even more advantages. The velocity of 
the p^y-Dirac cone is much larger than that of the s-band 
due to a much larger band width as shown in Fig. [5] The 
large energy scale here renders quantum degeneracy and 
the low temperature regime much more accessible. 

B. Strong coupling regime at n > | 

The physics in the strong coupling regime is very dif- 
ferent from that in the weak coupling regime. Much more 



crystalline ordered states appear in the strong coupling 
regime at commensurate fillings exhibiting rich structures 
of dimerization and trimerization orders. The relation of 
the filling (n) vs. at U/t = 10 is depicted in Fig. [T3] 
B. A scries of plateaus occur at commensurate fillings 
of (n) = | (i = 1 ~ 6), which correspond to a set of 
charge and bond ordered insulating states. The charge 
gap for each insulating state is at the order of ij| except 
for that at (n) = 1 which is at the order of U. Since 
these gapped state appears at strong interaction regime, 
they are not sensitive to a small t±. The band structure 
described in previous sections is completely changed by 
the strong interactions. Roughly speaking, at (n) > i, 
the preoccupied plaquette states exert strong effects to 
the extra particles and vice versa. The remaining part of 
the flat band disappear and the Dirac cone structure is 
also destroyed. 

At (n) = i, the strong coupling crystalline order- 
ing pattern is different from that in the weak coupling 
regime depicted in Fig. [T2J The system exhibits trimer- 
ized pattern as depicted in Fig. [TU Each trimer is re- 
pented as two thickened bonds and contains one par- 
ticle. In other words, each hexagon plaquette in the 
(n) = jr case are occupied by two particles. Such a 
state is also three-fold degenerate, and the other two 
equivalent states can be obtained by translating the sys- 
tem one lattice constant to right and left directions. Let 
us consider one plaquette unit cell and describe the or- 
bital configuration. We mark the six sites as 1 ~ 6. 
The p-orbital configurations at sites 1, 6 and 2 are p x , 
cos 8p x ± sin 9py(9 = 158.4°), respectively, and those at 
4, 5 and 3 are related by a reflection operation respect to 
x-axis. In other words, the occupied p-orbital at site 1 is 
parallel to x-axis, and that at site 6 is almost along the 
direction of bond (1,6) with a small deviation of 8.4°. 
The particle density at each site is ne,2, 3, 5 = 0.27 and 
ni,4 — 0.46. The bonding strength between neighouring 
sites i, j are defined as B^ = —((pi ■ &ij)(pj • e^) + h.c.) . 
There are four non-equivalent bond strengths of = 
(2, 1); (2,3); (2,5'); (1,4'), where 4' and 5' are the equiv- 
alent sites of 4 and 5 in the neighbouring plaquettes. 
We have B 2A = 0.58i y , B 2 ^ = 0.04i , , B 2 ,w = 0.14ty 
and B\ A ' = 0. The average bonding energy per site 
can be evaluated from the above bonding strengths as 
0.446in. Instead of the above trimer pattern, one might 
also think of the dimer covering with filling i in which 
only two third of sites are covered by dimers. However, a 
rough estimation of the average bonding energy per site 
is approximately about |iy, which is less energetically fa- 
vorable because particles are more localized in the dimer 
configuration. 

The crystalline order pattern at the filling of (n) = | 
is similar to that at (n) = i with each trimer containing 
two particles. In this case, the parameters above change 
to 9= 153°; n 6 . 2 , 3 .5 = 0.76, n 1A = 0.49; B 2 .i = OJOiy, 
B 2 ,3 = Q.07£||, B 2 ,5> = a.l2*|| and B lA , = 0.02fy. 

At (n) = i, the system exhibits a dimerized pattern 
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similar to that of (n) = g in the weak coupling regime as 
illustrated in Fig. [TO] The major difference is that the 
dimerized state here is an imcompressible insulating state 
while that in the weak coupling regime is with a small 
but still non-vanishing compressibility. The dimer is rep- 
resented by a thickened bond in which one particle hops 
back and forth. It can be considered as a superposition 
of the two states of two sites where one is occupied and 
the other is empty. There are only two non-equivalent 
bonding strength: £? lj6 = 0.95tj and B 12 = O.liii. The 
former is about one order larger that the latter, thus the 
system is in the strong dimerization limit. As shown in 
Fig. [12] B, the energy scale of this dimerized phase is 
set by f || , which is much larger than the usual one in 
dimerized magnetic systems with t 2 /U. 

The low energy physics in the dim er p hase should be 
described by a quantum dimer model [4l[ , which includes 
the quantum resonance of different patterns of dimer cov- 
erings. Although in the physical parameter regime the 
dimer crystal configuration in Fig. [13] is stabilized, it 
would be interesting to further investigate how to en- 
hance quantum fluctuations to achieve the quantum dis- 
ordered dimer liquid phase. The corresponding possible 
orbital liquid state in the p^-orbital systems would be 
an exciting state for a future study @ . 

The ordering pattern at another commensurate filling 
of (n) = | as shown in Fig. [TO] The p-orbital configura- 
tions at sites 1, 6 and 2 are p y , cos 9p x ± sin 9p y respec- 
tively with 9 = 150.2°, and those at 4, 5 and 3 are related 
by a reflection operation respect to £-axis. The particle 
density at each site is 716,2,3,5 = 0.94 and 711,4 = 0.62. 
The four non-equivalent bond strengths read B2.1 = 0.36, 
B 2 ,3 = 0.08, B 2 ,5< = 0.08 and B 1A > = 0.83. The |-filling 
state can be considered as doping the insulating state of 
one particle per site with i holes. Holes are mainly con- 
centrated on the positions of sites 1 and 4 in each unit 
cell. The corresponding bonds have the largest bonding 
strength. Such a state is the dimerized state of holes. 



For the close pack hexagon state at (n) = g, Eq. [35l 
can be easily calculated. We have (n(r) t ) — (^) 3 \i ; {k)\ 2 , 
where k = mrj (ht), and ip(k) is the Fourier transform of 
the plaquette-Wannier state depicted in Fig. [4] Thus 



N m -» -> 

C t (f,?) = T -(-)^( k )\ 2 \^(k')\ 2 



E 



S(k-k' - G'), 



(36) 



where '— ' ('+') is for fermions (bosons) respectively, G — 
mGi + nG' 2 with m, n integers, and k! = mf / (fit). After 
a spatial averaging and normalization, we find 



C t (d) = df 



C t (i 



d zf _ d 
2 > 2 



OC 



(n(f+l)) t (n(r-l))t 
T$>(fc-G'), 



(37) 



where k = rnd/(ht). All the <5-peaks are with equal 
weight because of the cancellation of the Fourier trans- 
form of the Wannier function, and the six-fold rotational 
symmetry in Fig. [4] 

For the crystalline ordering of fermions at other com- 
mensurate fillings, the noise correlation functions still ex- 
hibit the (5-peaks located at the same reciprocal wavevec- 
tors of the reduced Brillouin zone. However the form 
factors are more complicated. Generally, (n(r)) and 
Ct(f, r 1 ) can be calculated as 

(n(r) t ) cx J2^)Mk){n\pt(k)Pv(km), 

[IIS 

c t (fy) cx -\J2^(k)M^\pl(k)Pu(k')\n)\ 2 

x ^2s(k-k' -G'), (38) 



VI. TIME OF FLIGHT SPECTRA 

Noise correlation has become an important method 
to detect the ordering in cold atom systems in optical 
lattices [42|, [43|. In all the Mott-insulating states at 
commensurate fillings described in figures of [TO] [TO] HU 
and 1151 the enlarged unit cell contains six sites form- 
ing a plaquette. They should exhibit themselves in the 
noise correlation of the time of flight (TOF) signals. In 
the presence of the charge and bond orders, the recip- 
rocal wavevector of the reduced Brillouin zone becomes 



G[ 

2/ 



.0) 



l -G, 



G 2 and G' 2 



/ -27T 27T ^ 



<3V3a> ' 3 Wi n 3 W ^ ^2 V 3 ^ a 

jGi — \&2i where G\^ are the reciprocal wavevectors 
for the original Brillouin zone. The correlation function 
is defined as 



C t (r,7) = (n(f)n(f)) t - (n(r)) t (n(7)) u (35) 
where t is the flying time. 



where the Greek indices p, and v denote the Wannier 
functions for the 12 p XtV orbitals in one plaquette. In 
particular, due to the loss of the six-fold lattice rota- 
tional symmetry in [T4l and [TBI the noise spectra of Ct(d) 
should show the reduced two fold rotational symmetry. 
Fig. [TO] still keeps the six-fold rotational symmetry, but 
the weight of the (S-functions should not be the same as 
in Eq. EH 



VII. CONCLUSION AND DISCUSSION 

In summary, we have proposed the laboratory analog 
simulation of a new kind of artificial graphcnc, unavail- 
able in nature, where the p^^-orbitals are the key, unlike 
the real graphene made of the p z -orbital. This switching 
of orbitals, as shown in this work, lead to novel strong 
correlation physics which cannot be studied in the corre- 
sponding solid state graphene systems. 
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We have shown the band structure of p^-orbital 
honeycomb lattices contains both Dirac cones and flat 
bands. Particle interactions stabilize various incompress- 
ible Wigner crystal-like states at commensurate fillings. 
In particular, we have described the exact many body 
ground state at (n) — i, which exhibits close packed 
hexagon plaquette order. Various charge and bond order- 
ings appear in the strong coupling regime at higher com- 
mensurate fillings. These states exhibit their patterns in 
the noise correlation of time of flight experiments. Taking 
into account the recent exciting experimental realization 
of the p-orbital bosons [l8| and the fact that the hon- 
eycomb optical lattices were experimentally constructed 
quite some time ago (32|, the p^y-orbital counterpart of 
graphene may be achieved in the laboratory in the near 
future. 

Let us compare the Wigner crystal states in the p x ,y~ 
orbital systems with those in the electron gas systems. 
Quantum Monte-Carlo simulations show that the Wigner 
crystal state is stable in the very low density regime at 
r s > 39 in two dimensions, where r s is the ratio be- 
tween the average inter-electron distance and the Bohr 
radius 27]. The long range Coulomb interactions dom- 
inate over the kinetic energy when r s is large. In con- 
trast, even the shortest range repulsive interaction can 
stabilize the crystal state in the p^-orbital honeycomb 
lattice due to the suppression of the kinetic energy by 
the band flatness. Wigner crystal state also occurs in 
the fractional quantum Hall systems due to the suppres- 
sion of kinetic energy by the magnetic field [4J, |45| . At 
low filling factors, crystalline ordered states energetically 
win over the Laughlin liquid state. It is also interesting to 
note the difference between our system and the p 2 -orbital 
system of graphene, where the characteristic ratio be- 
tween Coulomb interaction and kinetic energy j^- (vf is 
the slope of the Dirac cone) is a constant independent of 



charge carrier density. As pointed out in Ref. [28,31,46], 
interactions in graphene are not strong enough to stabi- 
lize the Wigner crystal state at any density. 

Many interesting problems still remain open for fur- 
ther exploration, and we will leave them in future publi- 
cations. For example, for the spinful fermions with repul- 
sive interactions, it is natural to expect ferromagnetism 
due to the flat band structure. It would be interesting 
to study the competition between antiferromagnetic ex- 
change and flat band ferromagnetism. If interactions are 
attractive, the pairing problem and the corresponding 
BCS-BEC crossover in the flat band might prove inter- 
esting. On the other hand, if we load bosons into the 
flat band beyond the density of (n) = ^, the frustration 
effect due to the band flatness to the superfludity is a 
challenging problem. Most intriguing is the possibility 
of exotic incompressible states analogous to the Laughlin 
liquid in fractional quantum Hall effect. These cannot be 
captured within the mean-field approximation used here 
for n > 1/6. If one could devise appropriate variational 
liquid states projected into the flat band, these could be 
compared energetically with the Wigner crystals found 
here. Given the richness and surprises encountered in 
the fractional quantum Hall effect, flat band physics in 
optical lattices appears rife with possibility. 
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